arXiv: 1507.07772vl [math.NA] 28Jul2015 


High order numerical methods for networks of hyperbolic conservation laws 
coupled with ODEs and lumped parameter models 


Raul Borsche a ’*, Jochen Kall a 

a Erwin Schrodinger Strafie , TU Kaiserslautern, Building 48, 67663 Kaiserslautern, Germany 


Abstract 

In this paper we construct high order finite volume schemes on networks of hyperbolic conservation laws 
with coupling conditions involving ODEs. We consider two generalized Riemann solvers at the junction, 
one of Toro-Castro type and a solver of Harten, Enquist, Osher, Chakravarthy type. The ODE is treated 
with a Taylor method or an explicit Runge-Kutta scheme, respectively. Both resulting high order methods 
conserve quantities exactly if the conservation is part of the coupling conditions. Furthermore we present a 
technique to incorporate lumped parameter models, which arise from simplifying parts of a network. The 
high order convergence and the robust capturing of shocks is investigated numerically in several test cases. 

Keywords: ADER, Network, hyperbolic conservation law, WENO, generalized Riemann problem, 
Coupling, ODE, Runge-Kutta, Lumped Parameter Models 


1. Introduction 

Networks of hyperbolic PDEs arise from the modeling of many different problems, e.g. water and 
wastewater networks [Him gas pipelines traffic flow Eua, simulation of blood flow [0 HO] (TT] or 

cell migration [12] . The description of such networks is based on one dimensional conservation laws along 
the edges and suitable coupling conditions at the nodes. The simplest type of coupling uses a set of algebraic 
relations routing the flow between the arcs of the network. In many of the above applications further coupling 
conditions arise in which an ODE is located in the junction e.g. buffers storage tanks,manholes 010121 
or the heart M- A wide class of such coupling also occurs when so called lumped parameter models are 
applied to parts of the network [141 fT0[ lllj 0] . These models arise from simplification of the flow on the 
edges in regions where a coarser modeling can be afforded. 

In all these applications fast and accurate numerical methods are needed. For the one dimensional flow 
along the edges a huge variety of classical solvers for hyperbolic conservation laws is available P~5l Ig. 
Especially numerical methods of high order accuracy, e.g. WENO, ADER and DG schemes EaESHH], 
achieve remarkably accurate solutions relative to their computational costs. For the application on networks 
however mainly first order schemes have been developed. Recently a high order Riemann solver for purely 
algebraic coupling conditions was presented in m- 

In the present article we introduce two approaches to high order methods for vertices that can involve 
ODEs in addition to algebraic coupling conditions. The first method follows the Toro-Castro approach of 
[20] 19j. After solving a classical first order coupling problem linear coupling conditions for the temporal 
derivatives of the states at the junction are considered. The ODE is incorporated into this procedure by 
inserting the full Taylor-expansion of its solution into the coupling conditions. The second method adapts 
the Harten, Enquist, Osher, Chakravarthy approach [20, which solves a series of classical nonlinear Riemann 
problems. These problems are considered at the supporting points in time of an explicit Runge Kutta scheme, 
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Figure 1: Edge orientation convention 


which is used for the discretization of the ODE. In the context of this solver we investigate an efficient high 
order solver for lumped parameter models. Here we aim to exploit the underlying network structure for the 
numerical method. 

This paper is organized as follows: First we formulate the problem, specify the coupling conditions and 
define the generalized Riemann problem at a junction. In section [3] we recall the first order solver for such 
a coupling problem. The high order method of Toro-Castro type is presented in Section |4j the HEOC 
type solver in section [5j Based on these we describe a modification suited to networks including lumped 
parameter models. Finally we present convergence studies and several numerical examples in section [7] to 
show the quality of the schemes presented. 

2. Formulation of the problem and coupling conditions 

A Network J\f = (£ , V) consists of a set of edges £ = {Ei, ..., E^} which connect the vertices of the set 
V = {Vi,..., Vm}- On each edge E i: i = 1,..., h the quantities u l (x , t) £ are governed by a hyperbolic 
conservation law of the form 


d t u l + d x f\u l ) = 0 (1) 

with the flux function f l : —x R di , time t £ R + and location x £ [0, Lj], 

In the following we consider a single vertex V and assume all edges to be oriented outwards, as depicted 
in Figure ??. Starting from this setup, networks of arbitrary shape can be easily constructed by elementary 
transformations, e.g. m ■ To improve readability we drop the index j for all quantities at the junction, i.e. 
w = w 1 , and the spatial dependencies of the ids, which are evaluated at x = 0. 

At the vertex V we now assume a coupling of mixed algebraic-ODE type, i.e. 

§{v}(t),...,u n {t),w) =0 , 

w = F(u 1 (t),...,u n (t),w ), (2) 

u l (t) = w*(0, t) . 

The algebraic coupling conditions are given by the function $ : &" =1 x l 1 -> K c for n connected edges 
and the ODE is defined by the flux F : (^)' , =1 xl 1 R ! . Coupling conditions of that type arise from 
the modeling of e.g. storage components as manholes [3] or reservoirs [22], queues [231 HI], as well as from 
representing parts of complicated networks by so called lumped parameter models mm®- 

In order to keep the number of coupling conditions constant, the eigenvalues of the Jacobians 
A], l = 1,..., di need to be bounded away from zero for all states considered 

Ai < ... < y di , |A?| > e > 0 VI = 1,..., di . (3) 

To ensure that the correct number of coupling conditions is provided and that the problem is well posed 
[25] . we require: 

det (V u i>I>(ug,..., wo)i? +,1 |... ..., Ug, wo)R +,n ) 7 ^ 0 , (4) 
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where R +,l = [r d ._ c . +1 \... \r l d .\ denotes the matrix of all eigenvectors of V u i/ 1 which belong to positive 
eigenvalues. 

Ci = #{Aj|Aj->0} (5) 

denotes the number of positive eigenvalues on edge i and c = X^=i c i defines the total number of coupling 
conditions prescribed by 

3. The generalized Riemann problem at a junction 

One central building block for the construction of schemes in the ADER framework is the generalized 
Riemann problem. In order to develop similar high order methods for networks, a detailed understanding 
of the generalized Riemann problem at the junction is required. If additionally an ODE is located at the 
node, we aim to split the problem into two separate ones. On the PDE side we are looking for high order 
approximations to the Godunov states at the boundaries of the PDE domains, while simultaneously evolving 
the ODE in the vertex one time step. 

In the following we will discuss two variants to tackle this problem. The first is based on the classical 
ADER approach of Toro-Castro accompanied by a Taylor method for the ODE. The second one utilizes a 
Harten, Enquist, Osher, Chakravartlry solver for the PDE and a Runge-Kutta scheme for the ODE. 

Definition 1. Generalized Riemann problem of order k at a junction: 

Consider a algebraic-ODE type coupling Q of n edges governed by (jTJ). We call such a coupling situation 
with given initial state of the ODE wq and polynomial Riemann data u l (x,0) of order k a Generalized 
Riemann problem of order k at the junction 

... ,u n {t),w(t)) = 0 , 

w = F{v}(t),... ,u n {t),w(t)), 

U l (t ) = 0 

k —1 i 

u l (x,o) = J2 pl ifi 
1 =0 

w( 0) = wo . 

Analogously to the classical Riemann problem, the states at the left boundary of the coupled edges u l g (t) = 
lim T _).o+ u’(0, t + r) are called Godunov states. 

3.1. Solving the classical Riemann problem at the junction 

Before considering the generalized Riemann problem we investigate the classical Riemann problem at a 
junction with an ODE, i.e. the setup of definition [T] with k = 1 

$(V(f),. ..,u n {t),w(t)) = 0 , w = F(u 1 (t),.. .,u n (t),w(t)), 

u l (t) = u*(0, t) u\x 1 Q)=Po 

u>(0) = w 0 ■ 

In the case without ODE, the Godunov states are constant in time, i.e. u g (t) = u l g ( 0) Vf > 0. Due to 
the presence of the ODE the state in the junction can vary over time and thus also the Godunov states 
u l 0 (t) change. However for solving the generalized Riemann problem we are only interested in the states at 
t -* 0+. 

Since the solution of the ODE is continuous in time [25] we have limt_ > .o+ ^it) = Wq- Knowing the initial 
state of the ODE, the problem at t = 0+ reduces to a classical Riemann problem at a junction. This we 
can solve with the help of the so called Lax-Curves mm- 
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Solving such a classical Riemann problem at a junction is equivalent to finding a set of states that 
fulfill the algebraic part of the coupling conditions while being possible Godunov states accessible from the 
Riemann data on each edge. In order to be such an accessible state they have to lie on the concatenated 
Lax curves anchored in the right states u* = p l 0 [5B]: 

Li (£, ■ ■., e Ci , K) = L di _ Ci+1 (£, •) O ... o LI (4, <) , 

where the number of curves and free parameters c, is given by ©■ The operator o denotes the concatenation 
in the last variable, i.e for two functions g and h 


g{£ i 


; ; ') ° +1 > • ■ • ) £m; — 5(Cl! • ■ • : C(; MC/+1 1 ■ ■ • > Cm; 3z)) . 


Therefore we have to solve the equations 

^(Ll(ey r ),L 2 g (e,u 2 r ),...,L n g (C,K),w 0 ) =0 ( 6 ) 

for the unknowns £® = (Ci; • ■ ■ ;Cc )■ The local solvability of this system for states close to u r is assured by 
condition © ESI- Once the parameters £* are known, the Godunov states can be determined by evaluating 
the concatenated Lax-curves 

«*(o) = L*(r,t4). 

Thus the complete set of states at the junction at t = 0+ is given by {wj(0),..., u g ( 0), wq}. 


4. Generalized Riemann solver of Toro-Castro type 


In the ADER framework the Toro-Castro approach provides a procedure to construct a polynomial in 
time approximating the solution of the generalized Riemann problem at the considered interface. This is 
achieved by splitting up the problem into one classical non linear Riemann problem and k — 1 linearized 
Riemann problems for the temporal derivatives of the Godunov states. An extension of this procedure to 
junctions without ODEs was presented in [Ill- 

Let a generalized Riemann problem at an ODE junction be given as in definition |T| As a first step 
we solve the zeroth order classical Riemann problem as described in section |3.1| using the zero order data, 
i.e. u l r = Pq . Note that once the states at t = 0+ are known, we can directly evaluate the ODE w(0) = 
F(Ug( 0), ..., Ug(0), w(0)), which already provides some information about the development of w(t). 

In a second step we aim to compute the temporal derivatives of the involved states. As in the classical 
ADER framework, we obtain governing equations for the derivatives by differentiating the conservation laws 
with respect to t and obtain 

<9 t (<9 t V) + V P{u l g )d x {dtu l ) + ’sources’ = 0 k = 1 ,..., k max . (7) 

The term ’sources’ encompasses everything that only depends on derivatives of degree k and less. It can be 
dropped, since it does not immediately act on the Godunov states [T5]. These governing equations for the 
temporal derivatives are linear hyperbolic systems. Therefore the corresponding Lax curves are linear as 
well. The concatenated Lax curves to a Riemann problem with d^u l r as states on the right hand side have 
the short form 


U q k {C,u r ) = d, 


H 


+ R + '% 


k i 


R + ’ 1 = Di- 


a+1 1 


FcU 


( 8 ) 


where r® denotes the eigenvector corresponding to the j-th eigenvalue A® of \7f l (u l g ). Note that ([3| still 
holds as the Jacobian is the same as in the Riemann problem of order zero. 

In order to obtain the coupling conditions for the temporal derivatives, we differentiate $ with respect 
to time. The first order derivative of $ reads 




• - - - > w)d t u l g + V w ®(ul ,..., u n g , w)d t w = 0 , 


1 = 1 
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where all quantities are evaluated at t = 0 and u l g denotes the Godunov state at x = 0 on edge i. By 
inserting the ODE d t w = F(u g ,..., u g , w) we obtain 

V U $(U g , ...,Ug, W)d t Ug + . . . , Ug, W) F , . . . , Ug , W) = 0 

'Vu®(u 1 g ,...,u™,w)d t u + '&i(u 1 g ,...,u™,w) = 0 . (9) 

The function 'hi only depends on states at t = 0 and not on any temporal derivative. Thus we have a linear 
system governing the temporal derivatives of the Godunov states at the junction. Here and in the following 
we assume the coupling conditions $ to be sufficiently differentiable, otherwise we consider the usage of high 
order schemes not appropriate. 

For the derivatives of orders k > 2 additional terms arise, which we summarize as 

V u $(Mg,.. .,u g ,w)d k u g + \h k(u g ,d t u g ,... ,d k ~ 1 u g ,w,d t w,.. .,d k w) = 0 . 

Note that these lower order terms can not be dropped as in the classical ADER framework, i.e. equation 
0. since we can not expect that it acts delayed in any form. It is easy to see that dR only depends on 
lower order derivatives of u g but still on all k derivatives of w. This we can reduce by using the derivative 
of the ODE itself 


d k t w = dt" 1 F(u g (t) 
= Z(u 1 dtu 1 


,u g 

k -1 
t 


(t),w) 


,u n g ,d t u n g ,...,dt y 


g , w, d t w,... ,d k 1 w) 


( 10 ) 


Thus we end up with a linear equation for d k u g which only requires derivatives of lower order than k 


V u $( 


Ug, ■. .,u g ,w)d k u g + ^k(u g ,dtu g ,... ,<9 t u g ,w,d t w, ...,d£ 


•j) = 0 . 


( 11 ) 


Finally we can use this expression to obtain all the temporal derivatives of u g iteratively, by starting with 
the first order equation ([9]) and successively solving |TT] ) in increasing order of k. 

To solve each of these systems (111, we proceed as in the zeroth order case. The temporal derivatives 
of u g are governed by a linear conservation law Q and coupled by a set of linear coupling conditions ©• 
Thus we insert the concatenated linear Lax-curves Q with the free parameter 

n 

£ • • •, U n g , w) (R + '% + d k K) +*k = 0 




,u”,w)R+>% 






= 0 . 


By introducing the notations 

a i = V ... ,u™,w)R + ’ 

& = (£ (i 4 n ) T , 

this linear system can be written as 


A = (ai|a 2 | ... | a n ) , 
d k u r = (dfuj. d k u 2 r .. 


d?<) J 


+ w)dtu r + '& k =0 . 


( 12 ) 


The matrix A is exactly the one in 0 considered for the well-posedness of the coupling conditions. Since 
we have det A^O, we can solve for the unknowns 


& = A- 1 (-V u $(uJ, W)d k t u r - 4-fc) . 


In order to evaluate this expression, the temporal derivatives of the states within the edges d k u r are needed. 
Analogous to the classical ADER approach we start with a WENO reconstruction of the spatial initial data. 
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Since we are at the boundary of a domain, we have to use an one-sided reconstruction of type m to obtain 
spatial derivatives d^u r at the junction. These we can transform into temporal derivatives d^u r using the 
Caucliy-Kowalewski or Lax-Wendroff procedure da w\- 

Note that carrying out the Cauchy-Kowalewski procedure on the basis of the zeroth order Godunov state 
is not strictly necessary, but does not cause any additional computational costs either. It is however necessary 
to make the scheme revert to the classical Toro-Titarev ADER scheme for one on one coupling, as we showed 
in mi- Once the values of £& are determined, we can evaluate the Lax-curves ([8]) to obtain the temporal 
derivatives of the Godunov states d^u g . With these we can directly build a polynomial approximation of 
u g (t) at the interfaces of the junction. 

An approximation to the states of the ODE can now be constructed easily. Since for the computation of 
d^Ug we already determined all the derivatives of w using equation (10), we just combine these to a Taylor 
series approximating w{t). 

The scheme derived from the Toro-Castro type solver for the generalized Riemann problem at the junction 
with ODE can be summarized as follows: 


1. Obtain GRP data at the junction via one-sided polynomial reconstruction. 

2. Solve the zeroth order Riemann problem at the junction, as described in section [3.11 

3. Apply the Cauchy-Kowalewski procedure to obtain temporal polynomials as input data. 

4. Solve the generalized Riemann problem at the junction, as described above. 

5. Approximate the fluxes across cell interfaces at the junction using the temporal derivatives of the 
Godunov states. 

6. Update the state w in the junction by applying a Taylor scheme of order k. 

7. Fill the ghost cells at the junction if needed. 

8. Run a high order finite volume scheme e.g. ADER to compute the fluxes across interior cell interfaces. 

A detailed description how to fill the ghost cells, if needed, is given in PS- 

At this point we note that the above procedure requires, besides the Cauchy-Kowalevsky procedure, the 
derivatives up to order k of the ODE as well as those of the coupling conditions. 


5. Generalized Riemann solver at an ODE junction Harten, Enquist, Osher, Chakravarthy 

type 

A popular alternative to the Toro-Castro approach for solving generalized Riemann problems is the 
Harten, Enquist, Osher, Chakravarthy solver [201 . Instead of solving the linear governing equations for the 
derivatives, it uses several classical Riemann problems at different points in time to achieve a high order 
approximation of the flux. Here we present a HEOC type solver for junctions, which is accompanied by a 
Runge-Kutta scheme for the ODE. 

In the classical HEOC approach, the solution of a generalized Riemann problem is approximated by 
the solutions of a sequence of classical Riemann problems. First we fix a quadrature rule according to the 
desired order of the scheme. Then the classical Riemann problem at t = 0 is solved. Based on its solution the 
spatial data on each side of the interface can be flipped into the time domain using the Cauchy-Kowalevsky 
procedure. These polynomials can be evaluated at the supporting points of the given quadrature rule, such 
that a series of classical Riemann problems arises. Their solutions serve as approximation of the states at 
the interface at these points in time. This information can be inserted into the quadrature rule to obtain 
a high order approximation of the fluxes at the interface. In Figure [2] a schematic representation of this 
procedure is shown. 

In the following we adapt this approach to generalized Riemann problems at junctions with ODEs. Before 
considering the fluxes of the PDE, we start with the numerical method for the ODE. Here we choose an 
explicit Runge-Kutta scheme which is at least accurate of order k. Usually the coefficients of RK-schemes 
are given in form of a Butcher array 
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Figure 2: HEOC schematic. 


such that for an ODE w = F(t, w) the update formula reads 


W n +1 = y n + At E b i k i 1 
j =i 


ki = F(t n + dAt, w n + At ^ aijkj) . 

j =i 


(13) 


At this point we note that each RK-scheme of order k naturally provides a quadrature formula of order k 
with the supporting points CjAt, such that 


pAt « 

/ F(r)dr = At^6 i F(f n + c i Af) + 0(At fc+1 ) . 
- 70 l=i 


In the following we will use exactly these intermediate time levels ti = qAI to set up the HEOC coupling 
procedure. 

Before considering higher order terms we have to solve the classical Riemann problem at t = 0, 

,P'S, wo) = 0 • 


Following the procedure described in 3.1 we obtain the values of Godunov states at t = 0, i.e. wi(0),... u"(0). 


Using these values we can transform the spatial data from a one-sided polynomial reconstruction into 
temporal data via Cauchy-Kowalevsky procedure. Thus for each connected edge we obtain a temporal 
polynomial of the form 


h —1 

"'max 


4(*)= E 


Pi 


i,time 


1 =0 


t l 

T\ 


as input data for the generalized Riemann problem at the junction. With these values available, we now 
aim to solve the classical Riemann problems at the time levels ti 




(14) 


In order to apply the technique of section 3.1 we have to provide some approximation for the value of w(ti). 
This is naturally provided by the RK-scheme as in the second formula of (13) 


l -1 

wi = w 0 + A t^apiki . 

i—1 


(15) 
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As in the classical RK methods, the value wi is not necessarily an approximation of very high order, but 
chosen in such a way that in the final update the desired order is obtained. Note that for the evaluation of 
the stages ki = F(ti, Ug(ti), ..., Ug(ti), wi) i = 1,..., l — 1 values of Wi and u g (ti) are needed. But since we 
have chosen an explicit RK-scheme, only data from the previous l — 1 stages is used. 


Inserting (15) into (14) we obtain 




,Ug(tl),Wl) = 0 


and can solve this classical Riemann problem for the Godunov states at t[. 

Once all Riemann problems are solved successively, we can compose the solutions to determine the fluxes 
of the conservation laws and to update the ODE 


/i 


1 

2 


k 


Efc/WO) 

i=i 


k 

w(t + At) = w(t) + Atbiki . (16) 

i=l 


The coefficients bi in both formulas are those of the RK-scheme (131 and the stages ki have been already 
computed for the formula (15). The complete scheme can be summarized as 


1. Obtain GRP data at the junction via one-sided polynomial reconstruction. 

2. Solve the zeroth order Riemann problem at the junction, as described in section [3.1| 

3. Apply the Cauchy-Kowalewski procedure to obtain temporal polynomials as input data. 

4. Solve the classical Riemann problems at the times ti as described above. 

5. Approximate the fluxes across cell interfaces at the junction using the Godunov states at the time 
levels ti. 

6. Update the state w in the junction by applying the RK-scheme of order k. 

7. Run a high order finite volume scheme to compute the fluxes across interior cell interfaces. 


One advantage of this approach is that neither the derivatives of the coupling conditions nor those of the 
ODE are needed. Thus the only symbolic manipulation necessary is the CK procedure, which is required 
for ADERs scheme inside the domain anyway. Furthermore we can solve the ODE with some classical 
RK-scheme, which is helpful especially for complicated or large ODEs e.g. those that arise from lumped 
parameter models. 

The main disadvantage of this approach is the higher computational costs, since several nonlinear Rie¬ 
mann problems have to be solved instead of just one nonlinear and a couple of linear ones. Furthermore we 
do not have enough information at the interfaces to fill possible ghost cells at the computational boundary. 
Thus for the reconstruction in the interior of the domain but close to the boundary we use a reconstruction 
with variable stencil lengths m- 


5.1. Conservation of quantities 

In many applications and their corresponding models some quantities are conserved in the complete 
network, e.g. the total mass [24] H IS [TO] . This is not only established via the conservation laws on the 
edges, but also due to a careful choice of the coupling conditions and the ODE in the junction. 

Since the conservation is guaranteed to be exact in the interior of the edges by the numerical scheme, 
it is desirable that also the coupling procedure is conservative. The first order junction solver in section 


3.1 is conservative, since the same Godunov states are used in the coupling conditions as well as for the 
computation of the fluxes across the interfaces. If the ODE is updated by an explicit Euler scheme, its 
update relies on the same values as the coupling condition. Clearly at the junction the conservation is only 
guaranteed up to the precision of the numerical method used to solve the nonlinear system arising from the 
coupling conditions. 

In case of no ODE in the junction it has been proven in m that the Toro-Castro approach also conserves 
the selected quantities. This proof can be easily modified such that it fits the current setting. We just have 
to take care that the ODE is updated with exactly the same numerical values as those arising in the coupling 
procedure. Therefore it is mandatory to use a Taylor scheme for the ODE. 







For the HEOC approach at any intermediate time level ti a classical Riemann Problem at a junction is 
considered. If now some quantity is conserved in the underlying system, at each of these Riemann problems 
the fluxes of the resulting Godunov states and the flux of the ODE balance exactly. Since we choose the 
identical bis for the flux integration and the RK-scheme (161, this also holds for the final updates in the 
PDEs and the ODE. 


5.2. Source-terms 

In many applications the conservation law ([I]) is replaced by a balance law via introducing source terms. 
These usually do not affect the coupling procedure i^j- If we have a numerical scheme at hand that is 
capable to treat the source terms properly and include the sources into the Cauchy-Kowalevsky -procedure, 
all the above methods can be applied. Since the lower order terms in 0 are dropped, the sources do neither 
change the Lax curves nor the governing equations of the higher order derivatives. 


6. Lumped Parameter Models 

In a lot of real world applications the dimensions of the network exceed the affordable computational 
effort, e.g. capillaries in the circulatory system. At the same time a detailed description of the flow is only 
needed in certain areas of the network. Therefore it is often convenient to describe some parts of the network 
by simpler models. A wide class of such reduced models are the so called lumped parameter models, which 
are used to describe e.g. the human circulatory system m Gj or § as networks gj. 

In this section we explain a process to construct high order schemes for hybrid models of networks 
containing hyperbolic conservation laws and lumped parameter models. To have access to such a process 
in an algorithmic framework is especially of interest in the context of dynamical switching between highly 
resolved and reduced models j 3 ]. 

Consider one edge in the network of length L with a conservation law 

d t u(t, x) + d x ( f(u(t , a;))) = 0 . 

Following the approach proposed in El. a lumped parameter model is obtained by averaging over the whole 
spatial domain, i.e. 



( * ^ 



d t 

- / u(t,x)d x 
L Jo 

1 

+ L 

L )) - 0)) 


V-^ ) 


='-U‘ 


=>U’ = F(U, U l g ,U r g ) . 


Thus the averaged state U is governed by a simple ODE with the unknown Godunov states at the left 
and right boundary U l g and U g . If such an edge is connected to a node in the network, these values will 
be determined when solving the associated coupling conditions. The coupling procedure from section |3.1| 
remains unchanged by the averaging process, but the Lax curve have now to be anchored at the only 
accessible state U. This implies, if complete sections of the network are simplified to lumped parameter 
models in the above manner, that the PDEs are not coupled to a simple ODE but to a differential algebraic 
equation (DAE). The components of the lumped region can be summarized as follows 


w' = F(w , Ug) 


$ (u g ,w) = 0 
u g - L(w,£) = 0 


ODE part, originates from lumped edges and ODE 
parts already present in vertices, 

Algebraic constraints stem from the coupling condi¬ 
tions of the vertices, 

Lax curve condition to connect the Godunov states 
with the internal states. 


(18) 

(19) 
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Note that here w contains the averaged states of all lumped edges in a connected area and all possible ODE 
states in the junctions of this region. 

From a numerical point of view the algebraic constraints (191 do enforce the usage of an appropriate 
solver of DAEs in the junction, e.g. modified RK schemes PHI- In the following however we want to present 
a technique to construct a high order solver using the underlying network structure of the LPM. 

Since the lumped areas might be very large or vary in time, due to switching between the models, we base 
the following procedure on the HEOC approach. As in section [5] we first apply polynomial reconstruction 
to obtain spatial data on the PDE-edges. For each vertex in the lumped parts of the network we solve the 
following zeroth order coupling problem as © 


... ,u l g ,w 0 ) = 0 with 



edge is equipped with a PDE 
edge model is lumped 


This yields the Godunov states on the PDE edges at t = 0, which are used to flip the data into time by 
Cauchy-Kowalevsky procedure, providing polynomials u*’ tlme (t). We repeat this step for each stage l of the 
RK-scheme in the HEOC approach and for each vertex in the lumped network, i.e. we solve the zeroth order 
coupling problem at ti = CiAt, 


§{v} g ,...,u l g ,wi) = 0 with 


.i.time 


(ti) 


If the edge in question is a PDE 
if the edge in question is lumped 


The states U'j are known from the previous stage of the RK-scheme and the new ones are obtained by 

i i 

Ul+i = U 0 + At^2 ai+i,ikY , w 3 l+x = w 0 + At^2 ai+ijk™ , 

i i 

where and kf are the components of the intermediate values hi for the lumped edges or the ODEs in 
the junctions respectively. 

With this procedure we can keep the full network structure, such that we can easily select any part of this 
network and switch between the simplified and the more accurate model. The scheme can be summarized 
by the following steps 

1. Obtain GRP data at the junction via one-sided polynomial reconstruction. 

2. Solve the zeroth order Riemann problem at the junction, as described in section [3.11 

3. Apply the Cauchy-Kowalewski procedure to obtain temporal polynomials as input data. 

4. Compute each stage of the RK-scheme as described above. 

5. Approximate the fluxes across cell interfaces at the junction using the Godunov states at the time 
levels t;. 

6. Update the state w in the junction by applying the RK-scheme of order k. 

7. Run a high order finite volume scheme to compute the fluxes across interior cell interfaces. 


6.1. Source-terms in lumped parameter models 

Source terms can be treated in a straight forward way, since they do not affect the coupling procedure. 
Only when particular steady states should be preserved, the averaging process in has to be adapted 
accordingly. In case of a balance law we obtain 

d t u(t, x ) + d x ( f(u{t , a;))) = S(u(t, a;)) 

=>dtlJ + j [/([/;) - f(U l g )\ = S(u(t , *)) « S(U) 
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Figure 3: Problems arising from naive lumping. Left: Scheme can become very unbalanced. Right: Conservation of mass 
violated. 


For example in the case of the shallow water equations (21) a huge variety of well-balanced numerical schemes 
has been developed in order to incorporate the bottom elevation in a suitable manner e.g. [291 GEa E]. In 
this particular situation the approximation made to obtain S reads as 


— ghd x b ~ —gh 


rb(L)-b( 0) 


( 20 ) 


where h denotes the water level, b the bottom elevation and g the gravitational acceleration. 

As illustrated in figure [3] an independent treatment of flux an source can not lead to a method preserving 
steady states. Therefore we have to apply a reconstruction technique to the values U, which follows the 
bottom topography. This so called hydrostatic reconstruction [30] is a well known tool in this context. Thus 
whenever bottom elevation is considered we will apply the hydrostatic reconstruction on lumped edges in 
order to capture the ’lake at rest’ correctly. 

A further detail concerning this particular steady state is also depicted in figure [3| If an initial condition 
for the PDE model is provided on the network, the initial conditions of the lumped parameter model are 
obtained by the averaging process in ©• Since in ( |20| ) the bottom elevation is linearized, possible errors 
have to be compensated in the initial values of U. In the following we will just add the missing amount of 
water artificially to the initial states of the LPM model by modifying the reconstruction to 

1/ = — y u(x)dx + J b(x)dx — ~(b(L) + b( 0))^ . 


7. Numerical Examples 

In this section we investigate the above presented numerical methods in different test cases. As conser¬ 
vation law along the edges we choose the shallow water equations 


d t h + d x q = 0 

d t q + d x ( ^ ^ghA = 0 . 


( 21 ) 


h denotes the depth of the water and q is the discharge in x direction, g = 9.81 is the gravitational 
acceleration. 
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Figure 4: Split circle network. 


As coupling 
so called ’equal 


conditions in the junctions we consider two different sets of equations. The first one is the 
heights’ coupling, which reads for n connected edges 


n 

^2 Qi = 0 , 

i= 1 

hi — hi = 0 2 < i < n . 


( 22 ) 


The first equation states the conservation of mass at the junction. The remaining n — 1 equations force all 
heights at the junction to be at the same level. 

The second set of coupling conditions involve an ODE at the junction. Here we consider a storage tank 
or manhole at the coupling point, which is modeled by 


hm 

Qm 


Qm 

A m 

gA„ 


14 

2 g h\ 





(23) 


The two states are the vertical water level in the tank h m and the discharge Q m flowing into the volume. 
The constant A m is the horizontal cross sectional area of the storage tank. It is important that this model 
is always accompanied by the following set of coupling conditions 


^ ' Qi T Qm — 0 




— t4j -{-hi ) — 0 2 < i < n 

2 g hf 


The first equation again ensures the conservation of total mass in the coupled system. The following n — 1 
equations state the equality of the so called hydraulic heads or energy levels. As this name indicates these 
conditions provide the conservation of the total energy at the junction in case of smooth solutions 0 C 2 . 
The total energy in the coupled system is conserved due to equation (231 in the storage model [2]. Here we 
further note that ( [23| does not depend on the choice of the related edge since all hydraulic heads coincide. 

In the following convergence studies, we use a simple network consisting of three edges and two nodes, as 
shown in Figure [I] All edges are of the same length L = 25 and in both nodes a storage tank with A m = 1 
is placed. The initial data used for convergence studies needs not only to be smooth along the edges, but 
also has to satisfy the coupling conditions and its temporal derivatives up to the order of the schemes to be 
investigated. In the following we choose lake at rest like states at the junctions and a sufficiently smooth 
transition between their two water levels. This can be achieved by < 7 *(a;, 0) = 0 i = 1,2,3 and polynomial 
states of degree 15 for hi which are determined by the constraints 


hi(0, 0) = 2 , d*hi(0,0) = 0 


hi(L, 0)=3, d*hi(L, 0) = 0 , k = 1,... ,7 . 
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Figure 5: Initial Data and solution for the convergence studies [7.1| and |7.2| PDE states on the edges to the left, ODE states 
in the vertices on the right side. 


The initial states of the ODEs are also at rest, i.e. w 1 (0) = (2, 0) T and w 2 (0) = (3,0) T . In Figure[5]we show 
initial states as well as the reference solution for the PDEs on the left side, on the right hand side are the 
states of the ODEs in the nodes for t £ [0,2.4]. The computation of the reference solution was performed 
by a scheme of order 8 on a grid of 800 cells per edge. 

In all numerical examples the time step At is synchronized in the complete network according to a CFL 
number 0.95. The stability bounds of the ODEs are always less restrictive than those of the PDEs. 

7.1. Convergence study Toro-Castro 

The first convergence tests we perform for the Toro-Castro solver which is described in section [4] In 
Figure [6] the L 1 and the L°°-e rror at t = 2.4 are plotted against the reciprocal of the cell width Air. 
Furthermore the L 2 -errors of the involved ODEs are shown. Additionally for selected orders k the errors 
and rates of convergence are given in Table [l] All numerical solutions converge with the expected order. 

7.2. Convergence study Harten, Enquist, Osher, Chakravarthy 

For the HEOC-coupling procedure we repeat the same test for schemes up to order k = 6. The L 1 - and 
L°° errors for the PDEs and the L 2 -errors of the ODEs are shown in in Figure [7] The precise values and 
the corresponding convergence rates are presented in Table [2] 

The errors are within the same range as those of the Toro-Castro method and all solutions converge with 
the predicted order. 

7.3. Convergence study Lumped Parameter Model 

The final convergence test addresses lumped parameter models. Therefore we consider a network con¬ 
sisting of six edges and four vertices as depicted in Figure [8j In the two junctions V 2 and V 4 equal height 
coupling is applied, while in V\ and V 3 two manholes with A m = 1 are located. The lumping of section [6] is 
applied to the whole network except i.e. the green colored region. 

As initial conditions we choose water at rest qt(x,0) = 0 i = 1,...,6 with the constant water level 
hi(x, 0) = 0, i = 2,..., 6 , except for the first edge. There we take as data a polynomial of degree 14 such 
that 

M^,0) = 5.3, 


MM) = 5 , 
^( 0 , 0 ) = 0 , 
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hi(L, 0) = 5 , 

d%h 1 (L,0)=0 Vfc = l,...,6 
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Figure 6: Convergence si ngly |7.11 Splitcircle TT solver in L 1 and L°° norm for the edges, L 2 norm for the ODEs. 



k = 

2 

k = 

4 

k = 

5 

k = 

6 

k = 

8 

PDE L 1 norm 

N 

L 1 

o L i 

L 1 

Ov 

L 1 

Oli 

L 1 

o L i 

L 1 

Ov 

50 

5.51e-02 


9.32e-04 


4.34e-04 


2.37e-04 


1.02e-04 


100 

1.44e-02 

1.94 

5.71e-05 

4.03 

1.52e-05 

4.83 

3.81e-06 

5.96 

1.71e-06 

5.89 

200 

3.79e-03 

1.92 

2.95e-06 

4.28 

4.04e-07 

5.24 

5.77e-08 

6.04 

3.45e-09 

8.96 

400 

9.61e-04 

1.98 

1.86e-07 

3.98 

1.05e-08 

5.27 

1.46e-09 

5.31 

1.67e-ll 

7.69 


PDE L°° norm 


N 

L°° 


L°° 

O l ° o 

L°° 

Ol- 

L°° 


L°° 

Ol°° 

50 

1.06e-02 


2.89e-04 


2.77e-04 


1.38e-04 


3.85e-05 


100 

3.59e-03 

1.56 

5.13e-05 

2.50 

1.21e-05 

4.51 

3.43e-06 

5.34 

2.23e-06 

4.11 

200 

1.01e-03 

1.83 

2.06e-06 

4.64 

4.17e-07 

4.86 

7.81e-08 

5.45 

5.85e-09 

8.57 

400 

2.58e-04 

1.96 

1.43e-07 

3.85 

1.05e-08 

5.30 

1.83e-09 

5.42 

4.90e-ll 

6.90 


ODE L 2 norm 


N 

L 2 

O l 2 

L 2 

0 L 2 

L 2 

O l 2 

L 2 

0 L 2 

L 2 

O l 2 

33 

2.31e-03 


5.62e-05 


1.02e-05 


4.02e-06 


3.35e-07 


64 

6.22e-04 

1.98 

3.27e-06 

4.30 

4.53e-07 

4.70 

6.64e-08 

6.19 

1.63e-09 

8.05 

120 

1.60e-04 

2.16 

1.94e-07 

4.50 

1.67e-08 

5.25 

1.02e-09 

6.64 

7.02e-12 

8.66 

237 

4.01e-05 

2.04 

1.17e-08 

4.13 

5.53e-10 

5.01 

1.56e-ll 

6.14 

2.35e-14 

8.37 


Table 1: Convergence study [TT] Convergence rates splitcircle with TT solver. 
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Figure 7: Convergence study fT. 2| Splitcircle HEOC solver in L 1 and L°° 


norm for the edges, L 2 norm for the ODEs. 



k = 

2 

k = 

4 

k = 

5 

k = 

6 

PDE L L norm 

N 

Lr 

Oli 

L 1 

O l x 

L 1 

0 L r 

Lr 

o L i 

50 

4.39e-02 


6.80e-04 


1.38e-04 


1.94e-04 


100 

1.13e-02 

1.95 

5.21e-05 

3.71 

5.09e-06 

4.76 

3.48e-06 

5.80 

200 

3.03e-03 

1.90 

2.97e-06 

4.13 

1.48e-07 

5.10 

1.74e-08 

7.64 

400 

7.69e-04 

1.98 

2.02e-07 

3.88 

3.63e-09 

5.35 

3.18e-10 

5.78 


PDE L°° norm 


N 

L°° 

Olo c 

L°° 

o L °° 

L°° 

0 L oo 

L°° 


50 

1.05e-02 


3.49e-04 


6.87e-05 


1.40e-04 


100 

3.22e-03 

1.70 

4.59e-05 

2.93 

5.97e-06 

3.52 

9.02e-06 

3.96 

200 

8.73e-04 

1.88 

1.62e-06 

4.83 

1.77e-07 

5.08 

2.65e-08 

8.41 

400 

2.22e-04 

1.98 

1.20e-07 

3.76 

3.52e-09 

5.65 

1.50e-09 

4.14 


ODE L 2 norm 


N 

L 2 

Ol 2 

L 2 

0 L 2 

L z 

0 L 2 

L 2 

Ol 2 

33 

2.86e-03 


1.06e-05 


7.86e-06 


6.32e-06 


64 

8.16e-04 

1.89 

1.25e-06 

3.24 

2.79e-07 

5.04 

3.85e-07 

4.22 

120 

2.13e-04 

2.13 

2.90e-07 

2.32 

9.12e-09 

5.44 

1.23e-09 

9.14 

237 

5.40e-05 

2.02 

2.08e-08 

3.87 

2.63e-10 

5.21 

8.86e-12 

7.25 


Table 2: Convergence study |7.2| Convergence rates Splitcircle HEOC solver. 
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Figure 8 : Diamond network. The blue vertices V\ and V 3 are energy conserving manholes, i.e. ODE vertices, the rest algebraic 
equal height coupling vertices. The green area is lumped resulting in a network of one remaining Edge E\ and an LPM vertex 
with a state of 14 components. 
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Figure 9: Convergence study Diamond LPM p773| initial data and solutions: From top to bottom: PDE initial data on Ei, 
ODE components of height type and ODE components of impulse type. 
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k = 

2 

k = 

4 

k = 

5 

k = 

6 

PDE L 1 norm 

N 

L 1 

Ol 1 

L 1 

0 L i 

Lr 

0 L i 

L 1 

Oli 

50 

2.51e-02 


3.37e-05 


1.20e-05 


2.47e-06 


100 

6.56e-03 

1.94 

2.37e-06 

3.83 

2.06e-07 

5.87 

3.38e-08 

6.19 

200 

1.54e-03 

2.09 

1.47e-07 

4.00 

4.31e-09 

5.58 

3.79e-10 

6.48 

400 

3.79e-04 

2.02 

8.95e-09 

4.04 

1.06e-10 

5.34 

4.76e-12 

6.32 

800 

7.19e-05 

2.40 

5.48e-10 

4.03 

4.47e-12 

4.57 

3.75e-12 

0.34 

PDE L°° norm 

N 

L°° 

Ol- 

L°° 

Ol- 

L°° 

Ol- 

L°° 

Ol- 

50 

4.82e-03 


5.07e-06 


1.89e-06 


2.69e-07 


100 

1.32e-03 

1.87 

3.36e-07 

3.91 

3.09e-08 

5.93 

6.03e-09 

5.48 

200 

4.47e-04 

1.56 

2.25e-08 

3.90 

7.85e-10 

5.30 

7.67e-ll 

6.30 

400 

1.39e-04 

1.69 

1.34e-09 

4.08 

1.52e-ll 

5.69 

9.12e-13 

6.39 

800 

2.89e-05 

2.26 

8.13e-ll 

4.04 

8.79e-13 

4.12 

4.85e-13 

0.91 


ODE L 2 norm 


N 

L' z 

0 L 2 

L' z 

Ol 2 

L 1 

Ol ^ 

L' z 

0 L 2 

109 

1.27e-03 


3.45e-07 


7.49e-09 


3.56e-09 


213 

3.14e-04 

2.09 

2.50e-08 

3.92 

5.00e-10 

4.04 

4.10e-ll 

6.66 

425 

7.80e-05 

2.02 

1.62e-09 

3.96 

1.52e-ll 

5.06 

4.64e-13 

6.49 

845 

1.94e-05 

2.02 

1.03e-10 

4.02 

3.55e-13 

5.47 

3.24e-14 

3.87 

1686 

4.85e-06 

2.01 

6.45e-12 

4.00 

2.31e-14 

3.95 

3.07e-14 

0.08 


Table 3: Convergence study diamond LPM |7.3| convergence rates. 
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Figure 11: Modified split circle. 


hold. In Figure [9] initial data and solutions are shown. For edge E\ only the solution at t = 0 and t = 7 are 
plotted, whereas the states of all ODEs are shown on the full time interval [0, 7]. 

The error plots are shown in Figure |10[ with the corresponding data in Table [3J The solutions in all 
components involved converge with the designed order. 

7-4■ Capturing of shock waves 

In regions of smooth states the advantages of higher order methods are clearly indicated by the order of 
convergence. This does not hold for discontinuous solutions. In the following example we want to investigate 
the stability and accuracy of high order methods near shock waves. 

Therefore we consider a modified split circle network as depicted in Figure El In all the nodes the 
coupling conditions for a storage tank with A m = 1 are applied. As initial conditions we choose the constant 
water levels With the following initial data: 

h 3 {x, 0) = h 2 (x, 0) = h A {x, 0) = 5 , h 3 (x , 0) = 6 

and q l (x, 0) = 0 i = 1,..., 4. Using this setup instead of reusing the normal split circle with Riemann initial 
data ensures a clean Riemann problem at the junction without the interference of intermediary states. 

The evolution of the states in the ODE of vertex V\ is depicted in figure Figure [12] When zooming in, we 
can see that the 6 -tlr order scheme on the coarse grid of 50 cells is closer to the reference solution computed 
on a grid of 200 cells compared to the first order scheme. Even though the initial data is non smooth, which 
reduces the order of convergence to one, the solution benefits from the high order treatment. 

In Figure [13] the solutions along edges Ei and £ 2 are shown. We observe that the shock emerging from 
the vertex profits greatly in sharpness from a higher order scheme as well. At t = 4 the shock emerging 
from V 2 into E 4 after passing through E\ and E 2 is shown in Figure [l4j Again the solution of the 6 -th order 
scheme is much better than its first order counterpart. 


7.5. Shocks and lumped parameter models 

In this test we consider a larger network of 32 edges and 24 nodes. As shown in Figure [l5j the network is 
of a tree like structure, inspired by the human circulatory system e.g. [5]. The first four edges £ 1; £ 2 , £ 3 , £4 
and the last four £ 2 g, £ 30 , £ 31 , £32 have a length of L = 25, for all remaining edges we choose L = 2.5. In 


the nodes equal height coupling conditions (22) are used. In order to investigate the influence of parameter 


lumping on the solution we model the lower half of the network by an ODE as described in Section [6j 

As initial conditions we choose water at rest on all edges. On £1 we impose the following Riemann 
Problem 


h 1 ( 0 ,x) = 


3 x < 18.5 
2 else 
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q\ 0 ,x) = 0 , 









Figure 12: Schock investigation [7~4| Solution of the ODE in V\ over time. 
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Figure 13: Schock investigation [i~T| Solution of the PDE on E 1 /E 2 at t = 2. 
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Figure 14: Schock investigation |7.4| Solution of the PDE on E 4 at t = 4. 
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Figure 15: Split and join network with lumping of the green area. 


21 























Figure 16: Shocks and LPM [7^5] Height components of the LPM model over time compared to the averaged solution of the 
full PDE simulation. From top to bottom: Es, Eiq, E 24 and E 28 - 


whereas the remaining edges continue these constant states 

m 2 (0, x) = (3, Of , id(0, x) = (2, Of Vi > 3 . 

Due to the symmetric nature of this setting and without applying the lumping to the lower part of the 
network, the solution in some branches of the network coincide . In this case we have the following identities 


u 9 

,17 


X 10 

,18 



„,3 

.A 



U 

= u , 


„,5 


- f 7 

..8 

U = 

= U 

= U - 


„ 11 

12 

..13 

..14 

= U = 

U 

= U 

= U 

„ .19 

..20 

..21 

..22 

= U = 

U 

= U 

= U 

„,25 

-.26 

..27 

..28 

U = 

U 

= U 

= U 


-.29 

..30 



U 

= U 



X 15 

,23 


x 16 

,24 


Therefore it suffices to look at one edge of each group and we can directly compare the solutions of the 
lumped part of the network to its PDE counterparts. 

In Figure[l6]we show the heights H in the LPM vertex and the averaged heights on a corresponding edge 
located on the upper half of the network. Analogously the momentum components can be seen in Figure 

im 
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Figure 17: Shocks and LPM |7.5| Impulse components of the LPM model over time compared to the averaged solution of the 
full PDE simulation. From top to bottom: Eg, E\q, E 24 and E 2 8- 
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Figure 18: Shocks and LPM[7~5] Solution on Eog after the LPM vertex compared to the reference solution on E 3 q. 


Due to the very coarse spatial resolution the states of the LPM model can not resolve the incoming shock 
wave accurately. Despite this strong diffusion caused by the model, the ODE captures the general behavior 
of the flow. 

The correct capturing of the wave speed can be observed in Figure [18] Here the solution on the edges 
E 2 g and E 30 are shown. This provides a comparison between the shock wave, that has passed the lumped 
branch of the network and wave transported by the PDE model. Clearly the solution of the PDE model can 
resolve much finer structures. 

7.5.1. Sources in Lumped Parameter Models 

To demonstrate the necessity of hydrostatic reconstruction in LPM models, we simulate the following 
initial value problem consisting of shallow water equations on a split circle network with bottom elevation. E 2 
and the two vertices are turned into a LPM. We use bottom profiles of linear, polynomial and trigonometrical 
type: 

bi(x) = 0.3^ b 2 {x) 

25 

Hi(x) = 3 + l[i 25i ! 25] (x) H 2 (x) 

Pi{x) = H 1 (x) - bi(x) p 2 {x) 

The results of the simulations are shown in figure [19] and [20] The solutions on the PDE edges are shown at 
t = 0.3, the state of the ODE over the complete time interval. In the interior of the domain we can observe 
the asymptotically well balancedness of ADER schemes as reported in [2Dj, the slightly stronger oscillations 
at the vertices are caused by the higher sensitivity of the coupling conditions amplifying the unavoidable 
errors stemming from the one sided polynomial reconstruction. The big oscillations however are avoided as 
expected. 

8. Conclusion 

High order GRP solvers for ODE vertices and LPM models of Toro-Castro and Harten, Enquist, Osher, 
Chakravarthy type are introduced in this work. Extensive tests showed that they indeed exhibit the high 
order of convergence they were designed for. Numerical examples showed that these technique can indeed be 
used to build very accurate and stable numerical methods for networks of conservation laws including vertices 
with ODEs and lumped parameter models. Especially the solver of Harten, Enquist, Osher, Chakravarthy 
type is a vast improvement in terms of applicability over the Toro-Castro approach introduced in nn. 



= 3 

= H 2 (x) - b 2 {x) 


&3(*) = °. 3 (|+sm(^|)) 
H 3 (x) = 3 

p 3 (x) = H 3 (x) - b 3 {x) 
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Figure 19: Well balanced LPM [7lTT1 H on E\. 
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Figure 20: Well balanced LPM [7lTTl q on E\. 
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